/*
 * 	This file is a modification of phylip's drawtree.c. It is also part 
 *  of MVD_GUI.
 *
 *  MVD_GUI is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MVD_GUI is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with MVD_GUI.  If not, see <http://www.gnu.org/licenses/>.
 *  This version copyright 2010 Desmond Schmidt
 */
/* In addition to the above, the following notice also applies */
/* Version 3.6.  Copyright (c) 1986-2004 by the University of Washington and
  Written by Joseph Felsenstein and Christopher A. Meacham.  Additional code
  written by Sean Lamont, Andrew Keefe, Hisashi Horino, Akiko Fuseki, Doug
  Buxton and Michal Palczewski.  Permission is granted to copy, distribute,
  and modify this program provided that (1) this copyright message is
  not removed and (2) no fee is charged for this program. */

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <syslog.h>
#include <stdlib.h>
#include "phylip.h"
#include "draw.h"
#include "drawtree.h"
#define MAXITERATIONS 100
#define MINIMUMCHANGE   0.0001
#define GAP 0.5


void improveNodeAngle(node *pToNode, double medianDistance)
{
  	double forcePerpendicular, distance, 
    xDistance, yDistance, angleRotate, sinAngleRotate, cosAngleRotate, 
    norminalDistance, leftLimit, rightLimit, limitFactor;
  	node *pPivot;

  	/* Limit factor determinte how close the rotation can approach the 
     absolute limit before colliding with other branches */
  	limitFactor = (double)4 / (double)5;
  	pPivot = pToNode->back;

  	xDistance = nodep[pPivot->index-1]->xcoord - nodep[pToNode->index-1]->xcoord;
  	yDistance = nodep[pPivot->index-1]->ycoord - nodep[pToNode->index-1]->ycoord;
  	distance = sqrt( xDistance*xDistance + yDistance*yDistance  );

  	/* convert distance to absolute value and test if it is zero */
  	if ( fabs(distance) < epsilon) 
  	{
    	angleRotate = (double)0;
  	}
  	else 
  	{
    	leftRightLimits(pToNode, &leftLimit, &rightLimit);
    	norminalDistance = distance / medianDistance;
    	forcePerpendicular = forcePerpendicularOnNode(pPivot, pToNode,
        	medianDistance);
    	angleRotate = forcePerpendicular / norminalDistance;
    	/* Limiting the angle of rotation */
    	if ( angleRotate > 0 && angleRotate > limitFactor * leftLimit)
    	{
      		/* Left */
      		angleRotate = limitFactor * leftLimit;
    	}
    	else if ( -angleRotate > limitFactor * rightLimit )
    	/* angleRotate < 0 && */ 
    	{
      		/* Right */
      		angleRotate = - limitFactor * rightLimit;
    	}
  	}
  	angleRotate = (double).1 * angleRotate;

  	sinAngleRotate = sin(angleRotate);
  	cosAngleRotate = cos(angleRotate);

  	tilttrav(pToNode,  
    	&(nodep[pPivot->index - 1]->xcoord), 
        &(nodep[pPivot->index - 1]->ycoord), 
        &sinAngleRotate, &cosAngleRotate);

  	polarizeABranch(pToNode,
        &(nodep[pPivot->index - 1]->xcoord),
        &(nodep[pPivot->index - 1]->ycoord));
}
/** 
 * improvtrav for n-body. 
 * POPStack is the stack that is currently being used (popped); 
 * PUSHStack is the stack that is for the use of the next round 
 * (is pushed now) 
 */
void improvtravn(node *pStartingSubNode)
{
	stackElemType *pPUSHStackTop, *pPOPStackTop, *pTempStack;
  	node *pSubNode, *pBackStartNode, *pBackSubNode;
  	double medianDistance;
  	int noOfIteration;

  	/* Stack starts with no element on it */
  	pPUSHStackTop = NULL;
  	pPOPStackTop = NULL;

  	/* Get the median to relate force to angle proportionally. */
  	medianDistance = medianOfDistance(root, true);

  	/* Set max. number of iteration */
  	for ( noOfIteration = (long)0; noOfIteration < maxNumOfIter; noOfIteration++) 
	{
    	/* First, push all subNodes in the root node onto the stack-to-be-used
       	to kick up the process */
    	pSubNode = pStartingSubNode;
    	pushNodeToStack(&pPUSHStackTop, pSubNode);
    	while(pSubNode->next != pStartingSubNode)
    	{
      		pSubNode = pSubNode->next;
      		pushNodeToStack(&pPUSHStackTop, pSubNode);
    	} 
    	while (true)
    	{
      		/* Finishes with the current POPStack; swap the function of the stacks if 
               PUSHStack is not empty */
      		if (pPUSHStackTop == NULL)
      		{
              	/* Exit infinity loop here if empty. */
              	break;
      		}
      		else
      		{
              	/* swap */
              	pTempStack = pPUSHStackTop;
              	pPUSHStackTop = pPOPStackTop;
              	pPOPStackTop = pTempStack;
      		}
      
      		while (pPOPStackTop != NULL)
      		{
           		/* We always push the pivot subNode onto the stack!  That's when we
            	pop that pivot subNode, subNode.back is the node we apply the force
            	to (ToNode). Also, when we pop a pivot subNode, always push all 
            	pivot subNodes in the same ToNode onto the stack. */
              	popNodeFromStack(&pPOPStackTop, &pSubNode);
      
              	pBackStartNode = pSubNode->back;
              	if (pBackStartNode->tip == true)
              	{
                	/* tip indicates if a node is a leaf */
                	improveNodeAngle(pSubNode->back, medianDistance);
      
              	}
              	else
              	{
                	/* Push all subNodes in this pSubNode->back onto the
                	 * stack-to-be-used, after poping a pivot subNode. If
                	 * pSubNode->back is a leaf, no push on stack. */ 
                	pBackSubNode = pBackStartNode;
                	/* Do not push this pBackStartNode onto the stack!  Or the
                	 * process will never stop. */
                	while(pBackSubNode->next != pBackStartNode)
                	{
                  		pBackSubNode = pBackSubNode->next;
                  		pushNodeToStack(&pPOPStackTop, pBackSubNode);
                	} 
                	/* improve the node even if it is not a leaf */
                	improveNodeAngle(pSubNode->back, medianDistance);
              	}
      		}
    	}
  	}
}
/* traverse tree trying different tiltings at each node */
void improvtrav(node *p)
{
  	double xx, yy, cosphi, sinphi;
  	double langle, rangle, sumrot, olddiff;
  	node *pp, *qq, *ppp;;

  	if (p->tip)
    	return;
  	xx = p->xcoord;
  	yy = p->ycoord;
  	pp = p->next;
  	do 
	{
    	leftrightangle(pp, xx, yy);
    	pp = pp->next;
  	} while ((pp != p->next));
  	if (p == root) 
	{
    	pp = p->next;
    	do 
		{
      		qq = pp;
      		pp = pp->next;
    	} while (pp != root);
    	p->righttheta = qq->righttheta;
    	p->lefttheta = p->next->lefttheta;
  	}
  	qq = p;
  	pp = p->next;
  	ppp = p->next->next;
  	do 
	{
    	langle = qq->righttheta - pp->lefttheta;
    	rangle = pp->righttheta - ppp->lefttheta; 
    	while (langle > pi)
      		langle -= 2*pi;
    	while (langle < -pi)
      		langle += 2*pi;
    	while (rangle > pi)
      		rangle -= 2*pi;
    	while (rangle < -pi)
      		rangle += 2*pi;
    	olddiff = fabs(langle-rangle);
    	sumrot = (langle - rangle) /2.0;
    	if (sumrot > langle)
      		sumrot = langle;
    	if (sumrot < -rangle)
      		sumrot = -rangle;
    	cosphi = cos(sumrot);
    	sinphi = sin(sumrot);
    	if (p != root) 
		{
      		if (fabs(sumrot) > maxchange)
        		maxchange = fabs(sumrot);
      		pp->back->oldtheta += sumrot;
      		tilttrav(pp->back,&xx,&yy,&sinphi,&cosphi);
      		polarize(pp->back,&xx,&yy);
      		leftrightangle(pp, xx, yy);
      		langle = qq->righttheta - pp->lefttheta;
      		rangle = pp->righttheta - ppp->lefttheta; 
      		while (langle > pi)
        		langle -= 2*pi;
      		while (langle < -pi)
        		langle += 2*pi;
      		while (rangle > pi)
        		rangle -= 2*pi;
      		while (rangle < -pi)
        		rangle += 2*pi;
      		while ((fabs(langle-rangle) > olddiff) && (fabs(sumrot) > 0.01)) 
			{
        		sumrot = sumrot /2.0;
        		cosphi = cos(-sumrot);
        		sinphi = sin(-sumrot);
        		pp->back->oldtheta -= sumrot;
        		tilttrav(pp->back,&xx,&yy,&sinphi,&cosphi);
        		polarize(pp->back,&xx,&yy);
        		leftrightangle(pp, xx, yy);
        		langle = qq->righttheta - pp->lefttheta;
        		rangle = pp->righttheta - ppp->lefttheta; 
        		if (langle > pi)
          			langle -= 2*pi;
        		if (langle < -pi)
          			langle += 2*pi;
        		if (rangle > pi)
          			rangle -= 2*pi;
        		if (rangle < -pi)
          			rangle += 2*pi;
      		}
    	}
    	qq = pp;
    	pp = pp->next;
    	ppp = ppp->next;
  	} while (((p == root) && (pp != p->next)) || ((p != root) && (pp != p)));
  	pp = p->next;
  	do 
	{
    	improvtrav(pp->back);
    	pp = pp->next;
  	}
  	while (((p == root) && (pp != p->next)) || ((p != root) && (pp != p)));
}
/* get width and depth beyond each node */
void getwidth(node *p)
{
  	double nw, nd;
  	node *pp, *qq;

  	nd = 0.0;
  	if (p->tip)
    	nw = 1.0;
  	else 
	{
	    nw = 0.0;
    	qq = p;
    	pp = p->next;
    	do 
		{
      		getwidth(pp->back);
      		nw += pp->back->width;
      		if (pp->back->depth > nd)
        		nd = pp->back->depth;
      		pp = pp->next;
    	} while (((p != root) && (pp != qq)) || ((p == root) && (pp != p->next)));
  	}
  	p->depth = nd + p->length;
  	p->width = nw;
}
/* compute x and y coordinates */
void coordtrav(node *p, double *xx, double *yy)
{
  	node *pp;

  	if (!p->tip) 
	{
    	pp = p->next;
    	while (pp != p) 
		{
      		coordtrav(pp->back, xx, yy);
      		pp = pp->next;
      		if (p == root)
        		coordtrav(p->back, xx, yy);
    	}
  	}
  	(*xx) = p->r * cos(p->theta);
  	(*yy) = p->r * sin(p->theta);
  	if ((*xx) > maxx)
    	maxx = (*xx);
  	if ((*xx) < minx)
    	minx = (*xx);
  	if ((*yy) > maxy)
    	maxy = (*yy);
  	if ((*yy) < miny)
    	miny = (*yy);
  	p->xcoord = (*xx);
  	p->ycoord = (*yy);
}
/* use angles calculation to improve node coordinate placement */
void coordimprov(double *xx, double *yy)
{
  	int i;

  	if (nbody)       /* n-body algorithm */
    	improvtravn(root);
  	else 
	{          /* equal-daylight algorithm */
	    i = 0;
	    do 
		{
      		maxchange = 0.0;
      		improvtrav(root);
      		i++;
    	} while ((i < MAXITERATIONS) && (maxchange > MINIMUMCHANGE));
  	}
}
/* compute coordinates for tree */
void calculate()
{
  	double xx, yy;
  	int i;
  	double nttot, fontheight, labangle=0, top, bot, rig, lef;

  	for (i = 0; i < nextnode; i++)
    	nodep[i]->width = 1.0;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->xcoord = 0.0;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->ycoord = 0.0;
  	if (!uselengths) 
	{
    	for (i = 0; i < nextnode; i++)
      		nodep[i]->length = 1.0;
  	} 
	else 
	{
    	for (i = 0; i < nextnode; i++)
      	nodep[i]->length = fabs(nodep[i]->oldlen);
 	}
  	getwidth(root);
  	nttot = root->width;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->width = nodep[i]->width * spp / nttot;
  	if (!improve)
    	plrtrans(root, treeangle, treeangle - ark / 2.0, treeangle + ark / 2.0);
  	else plrtrans(root, treeangle, treeangle - pi, treeangle + pi);
  		maxx = 0.0;
  	minx = 0.0;
  	maxy = 0.0;
  	miny = 0.0;
  	coordtrav(root, &xx,&yy);
  	fontheight = heighttext(fontname);
  	if (labeldirec == fixed)
    	labangle = pi * labelrotation / 180.0;
  	textlength = (double*) Malloc(nextnode*sizeof(double));
  	firstlet = (double*) Malloc(nextnode*sizeof(double));
  	for (i = 0; i < nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		textlength[i] = lengthtext(nodep[i]->nayme, nodep[i]->naymlength,
				fontname);
      		textlength[i] /= fontheight;
      		firstlet[i] = lengthtext(nodep[i]->nayme,1L,fontname)
                      / fontheight;
    	}
  	}
  	if (spp > 1)
    	labelheight = charht * (maxx - minx) / (spp - 1);
  	else
    	labelheight = charht * (maxx - minx);
  	if (improve) 
	{
		coordimprov(&xx,&yy);
		maxx = 0.0;
		minx = 0.0;
		maxy = 0.0;
		miny = 0.0;
		coordtrav(root, &xx,&yy);
  	}
  	topoflabels = 0.0;
  	bottomoflabels = 0.0;
  	rightoflabels = 0.0;
  	leftoflabels = 0.0;
  	for (i = 0; i < nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		if (labeldirec == radial)
        		labangle = nodep[i]->theta;
      		else if (labeldirec == along)
        		labangle = nodep[i]->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		if (cos(labangle) < 0.0 && labeldirec != fixed)
        		labangle -= pi;
      		firstlet[i] = lengthtext(nodep[i]->nayme,1L,fontname)
				/ fontheight;
      		top = (nodep[i]->ycoord - maxy) / labelheight + sin(nodep[i]->oldtheta);
      		rig = (nodep[i]->xcoord - maxx) / labelheight + cos(nodep[i]->oldtheta);
      		bot = (miny - nodep[i]->ycoord) / labelheight - sin(nodep[i]->oldtheta);
      		lef = (minx - nodep[i]->xcoord) / labelheight - cos(nodep[i]->oldtheta);
      		if (cos(labangle) * cos(nodep[i]->oldtheta) +
          		sin(labangle) * sin(nodep[i]->oldtheta) > 0.0) 
			{
        		if (sin(labangle) > 0.0)
          			top += sin(labangle) * textlength[i];
        		top += sin(labangle - 1.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) < 0.0)
          			bot -= sin(labangle) * textlength[i];
        		bot -= sin(labangle - 0.75 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			rig += cos(labangle - 0.75 * pi) * GAP * firstlet[i];
        		else
          			rig += cos(labangle - 1.25 * pi) * GAP * firstlet[i];
        		rig += cos(labangle) * textlength[i];
        		if (sin(labangle) > 0.0)
          			lef -= cos(labangle - 1.25 * pi) * GAP * firstlet[i];
        		else
          			lef -= cos(labangle - 0.75 * pi) * GAP * firstlet[i];
      		} 
			else 
			{
        		if (sin(labangle) < 0.0)
          			top -= sin(labangle) * textlength[i];
        		top += sin(labangle + 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			bot += sin(labangle) * textlength[i];
        		bot -= sin(labangle - 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			rig += cos(labangle - 0.25 * pi) * GAP * firstlet[i];
        		else
          			rig += cos(labangle + 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) < 0.0)
          			rig += cos(labangle) * textlength[i];
        		if (sin(labangle) > 0.0)
          			lef -= cos(labangle + 0.25 * pi) * GAP * firstlet[i];
        		else
          			lef -= cos(labangle - 0.25 * pi) * GAP * firstlet[i];
        		lef += cos(labangle) * textlength[i];
      		}
      		if (top > topoflabels)
        		topoflabels = top;
      		if (bot > bottomoflabels)
        		bottomoflabels = bot;
      		if (rig > rightoflabels)
        		rightoflabels = rig;
      		if (lef > leftoflabels)
        		leftoflabels = lef;
    	}
  	}
  	topoflabels *= labelheight;
  	bottomoflabels *= labelheight;
  	leftoflabels *= labelheight;
  	rightoflabels *= labelheight;
}
/* compute coordinates of tree for plot or preview device */
void rescale()
{
  long i;
  double treeheight, treewidth, extrax, extray, temp;

  treeheight = maxy - miny + topoflabels + bottomoflabels;
  treewidth = maxx - minx + rightoflabels + leftoflabels;
  if (grows == vertical) {
    if (!rescaled)
      expand = bscale;
    else {
      expand = (xsize - 2 * xmargin) / treewidth;
      if ((ysize - 2 * ymargin) / treeheight < expand)
        expand = (ysize - 2 * ymargin) / treeheight;
    }
    extrax = (xsize - 2 * xmargin - treewidth * expand) / 2.0;
    extray = (ysize - 2 * ymargin - treeheight * expand) / 2.0;
  } else {
    if (!rescaled)
      expand = bscale;
    else {
      expand = (ysize - 2 * ymargin) / treewidth;
      if ((xsize - 2 * xmargin) / treeheight < expand)
        expand = (xsize - 2 * xmargin) / treeheight;
    }
    extrax = (xsize - 2 * xmargin - treeheight * expand) / 2.0;
    extray = (ysize - 2 * ymargin - treewidth * expand) / 2.0;
  }
  for (i = 0; i < (nextnode); i++) {
    nodep[i]->xcoord = expand * (nodep[i]->xcoord - minx + leftoflabels);
    nodep[i]->ycoord = expand * (nodep[i]->ycoord - miny + bottomoflabels);
    if (grows == horizontal) {
      temp = nodep[i]->ycoord;
      nodep[i]->ycoord = expand * treewidth - nodep[i]->xcoord;
      nodep[i]->xcoord = temp;
    }
    nodep[i]->xcoord += xmargin + extrax;
    nodep[i]->ycoord += ymargin + extray;
  }
}
void plotlabels( char *fontname )
{
	int i;
	double compr, dx = 0, dy = 0, labangle, sino, coso, cosl, sinl,
		cosv, sinv, vec;
  	boolean right;
  	node *lp;

  	compr = xunitspercm / yunitspercm;
  	if ( penchange == yes )
    	changepen(labelpen);
  	for (i = 0; i<nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		lp = nodep[i];
      		labangle = labelrotation * pi / 180.0;
      		if (labeldirec == radial)
        		labangle = nodep[i]->theta;
      		else if (labeldirec == along)
        		labangle = nodep[i]->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		if (cos(labangle) < 0.0)
        		labangle -= pi;
			sino = sin(nodep[i]->oldtheta);
			coso = cos(nodep[i]->oldtheta);
			cosl = cos(labangle);
			sinl = sin(labangle);
      		right = ((coso*cosl+sino*sinl) > 0.0) || (labeldirec == middle);
      		vec = sqrt(1.0+firstlet[i]*firstlet[i]);
      		cosv = firstlet[i]/vec;
      		sinv = 1.0/vec;
      		if (labeldirec == middle) 
			{
          		if ((textlength[i]+1.0)*fabs(tan(nodep[i]->oldtheta)) > 2.0) 
				{
            		dx = -0.5 * textlength[i] * labelheight * expand;
            		if (sino > 0.0) 
					{
              			dy = 0.5 * labelheight * expand;
              			if (fabs(nodep[i]->oldtheta - pi/2.0) > 1000.0)
                			dx += labelheight * expand / (2.0*tan(nodep[i]->oldtheta));
            		} 
					else 
					{
              			dy = -1.5 * labelheight * expand;
              			if (fabs(nodep[i]->oldtheta - pi/2.0) > 1000.0)
                			dx += labelheight * expand / (2.0*tan(nodep[i]->oldtheta));
            		}
          		}
          		else 
				{
            		if (coso > 0.0) 
					{
              			dx = 0.5 * labelheight * expand;
              			dy = (-0.5 + (0.5*textlength[i]+0.5)*tan(nodep[i]->oldtheta))
                    		* labelheight * expand;
            		}
            		else 
					{
              			dx = -(textlength[i]+0.5) * labelheight * expand;
              			dy = (-0.5 - (0.5*textlength[i]+0.5)*tan(nodep[i]->oldtheta))
                    		* labelheight * expand;
            		}
          		}
      		} 
			else 
			{
        		if (right) 
				{
					dx = labelheight * expand * coso;
					dy = labelheight * expand * sino;
					dx += labelheight * expand * 0.5 * vec * (-cosl*cosv+sinl*sinv);
					dy += labelheight * expand * 0.5 * vec * (-sinl*cosv-cosl*sinv);
        		} 
				else 
				{
					dx = labelheight * expand * coso;
					dy = labelheight * expand * sino;
					dx += labelheight * expand * 0.5 * vec * (cosl*cosv+sinl*sinv);
					dy += labelheight * expand * 0.5 * vec * (sinl*cosv-cosl*sinv);
					dx -= textlength[i] * labelheight * expand * cosl;
					dy -= textlength[i] * labelheight * expand * sinl;
        		}
      		}
      		plottext( lp->nayme, lp->naymlength,
				labelheight * expand * xscale / compr, compr,
				xscale * (lp->xcoord + dx + xoffset),
				yscale * (lp->ycoord + dy + yoffset), -180 * labangle / pi,
				fontname );
		}
	}
	if (penchange == yes)
		changepen(treepen);
} 
/**
 * Plot part or all of tree on the plotting device 
 */
void plottree( node *p, node *q )
{
	double x1, y1, x2, y2;
	node *pp;

	x2 = xscale * (xoffset + p->xcoord);
	y2 = yscale * (yoffset + p->ycoord);
	if (p != root) 
	{
		x1 = xscale * (xoffset + q->xcoord);
		y1 = yscale * (yoffset + q->ycoord);
		plot(penup, x1, y1);
		plot(pendown, x2, y2);
  	}
  	if (p->tip)
    	return;
  	pp = p->next;
  	do 
	{
		plottree(pp->back, p);
		pp = pp->next;
  	} while (((p == root) && (pp != p->next)) 
		|| ((p != root) && (pp != p)));
}
/** Initialize parameters */
void initialparms()
{
	paperx = 20.04;
	pagex  = 20.04;
	papery = 28.74;
	pagey  = 28.74;
	plotrparms(spp);
	grows = vertical;
	treeangle = pi / 2.0;
	ark = 2 * pi;
	regular = false;
	rescaled = true;
	bscale = 1.0;
	labeldirec = middle;
	xmargin = 0.08 * xsize;
	ymargin = 0.08 * ysize;
	labelrotation = 0.0;
	plotter = DEFPLOTTER;
	hpmargin = 0.02*pagex;
	vpmargin = 0.02*pagey;
	labelavoid = false;
}
/**
 * Initializes a node
 */
void initdrawtreenode( node **p, node **grbg, node *q, int len,
 	int nodei, int *ntips, int *parens, initops whichinit,
	pointarray treenode, pointarray nodep, char *str, char *ch,
	char *inTree, int *i )
{
	int j;
  	boolean minusread;
  	double valyew, divisor;

  	switch ( whichinit ) 
	{
	  	case bottom:
			gnu(grbg, p);
			(*p)->index = nodei;
			(*p)->tip = false;
			for (j=0;j<MAXNCH;j++)
			  (*p)->nayme[j] = '\0';
			nodep[(*p)->index - 1] = (*p);
			break;
	  	case nonbottom:
			gnu(grbg, p);
			(*p)->index = nodei;
			break;
	  	case tip:
			(*ntips)++;
			gnu(grbg, p);
			nodep[(*ntips) - 1] = *p;
			setupnode(*p, *ntips);
			(*p)->tip        = true;
			(*p)->naymlength = len ;
			strncpy ((*p)->nayme, str, MAXNCH);
			break;
	  	case length:
			processlength(&valyew, &divisor, ch,
			  &minusread, inTree, i, parens);
			if (!minusread)
			  (*p)->oldlen = valyew / divisor;
			else
			  (*p)->oldlen = fabs(valyew/divisor);
			if ((*p)->oldlen < epsilon)
			  (*p)->oldlen = epsilon;
			if ((*p)->back != NULL)
			  (*p)->back->oldlen = (*p)->oldlen;
			break;
	  	case hsnolength:
			haslengths = false;
			break;
	  	default:        /* cases hslength,iter,treewt,unitrwt        */
			break;        /* should not occur                        */
  	}
}
/**
 * Prepare for drawing the tree
 */
static void setUpEnvironment()
{
	node *q, *r;
  	boolean firsttree;

	grbg = NULL;
  	pointarray treenode = NULL;   
  	if ( strlen(inTreePath) > 0 )
	{
		intree = fopen( inTreePath, "r" );
	}
	else
	{
		intree = stdin;
	}
  	if ( intree != NULL )
	{
		char *treeContents;
		readTreeContents( intree, &treeContents );
		firsttree = true;
		allocate_nodep( &nodep, treeContents, &spp );
  		treeread( treeContents, &root, treenode, &goteof, &firsttree, nodep, &nextnode,
      		&haslengths, &grbg, initdrawtreenode, true, -1 );
		if ( !haslengths )
			uselengths = false;
		if ( intree != stdin )
			fclose( intree );
		q = root;
  		r = root;
		while (!(q->next == root))
			q = q->next;
  		q->next = root->next;
  		root = q;
  		chuck(&grbg, r);
  		nodep[spp] = q;
  		where = root;
  		rotate = true;
  		initialparms();
	}
	else
		syslog(LOG_ERR,"Couldn't open file %s\n",inTreePath );
}

/**
 * Parse the argument list and set globals
 * Currently: -i for intree (path), -p for plotfile (path)
 * and -n for max iterations (integer), -f for a postscript 
 * font-name No arguments are required. If unspecified, the 
 * infile is stdin and the outfile is stdout. Default font
 * is Times-Roman.
 * @return 1 if the args were OK, 0 otherwise
 */
static int readArgs( int argc, char ** argv )
{
	int i;
	int temp;
	for ( i=1;i<argc;i++ )
	{
		if ( strlen(argv[i]) > 1 )
		{
			if ( argv[i][0] == '-' && argc > i+1 )
			{
				switch ( argv[i][1] )
				{
					case 'i':
						strncpy(inTreePath,argv[i+1],FILE_NAME_LEN);
						break;
					case 'p':
						strncpy(plotFilePath,argv[i+1],FILE_NAME_LEN);
						break;
					case 'n':
						sscanf( argv[i+1], "%i", &maxNumOfIter );
						if ( maxNumOfIter > 0 )
							maxNumOfIter = 50;
						break;
					case 'f':
						strncpy( fontname, argv[i+1], FONT_NAME_LEN );
						break;
					case 'B':
						if (strcmp(argv[i+1],"true")==0)
      						uselengths = true;
    					else
							uselengths = false;
						break;
					case 'C':
						sscanf(argv[i+1],"%lf%*[^\n]", &charht);
						break;
					case 'I':
						sscanf(argv[i+1],"%d", &temp);
						// no improvement
						if ( temp == 0 )
						{
							nbody = false;
							improve = false;
						}
						// equal daylight
						else if ( temp == 1 )
						{
							nbody = false;
							improve = true;
						}
						// nbody
						else if ( temp == 2 )
						{
							nbody = true;
							improve = true;
						}
						break;
					default:
						return 0;
						break;
				}
			}
		}
	}
	return 1;
}
/**
 * Set sensible defaults for everything
 */
static void setDefaults()
{
	inTreePath[0] = 0;
	plotFilePath[0] = 0;
	maxNumOfIter = 50;
	strncpy( fontname, "Times-Roman", FONT_NAME_LEN );
	uselengths = haslengths;
	charht = 0.3333;
	improve = true;
	nbody = false;
}
/** 
 * Originally: loop to make preview window and decide what to do with it 
 * Intially cut out but then realised that it made initialisations.
 */
void user_loop()
{
    xscale = xunitspercm;
    yscale = yunitspercm;
    plotrparms(spp);
    calculate();
    rescale();
}
/**
 * Main entry point. Completely rewritten because we don't want 
 * (a) dependencies on non-ANSI-C software such as X11, and (b) 
 * to remove the user interaction. Just a plain old commandline 
 * interface with options works best for integration with an 
 * external GUI. Also we will read from stdin and write to stdout 
 * so we can pipe this command in and out of others such as 
 * ghostscript without writing any files.
 */
int main( int argc, char **argv )
{
	setDefaults();
	if ( readArgs(argc,argv) )
	{
		setUpEnvironment();
		user_loop();
		if ( strlen(plotFilePath)>0 )
		{
			plotfile = fopen( plotFilePath, "w" );
			if ( plotfile == NULL )
			{
				syslog(LOG_ERR,"Couldn't open %s for writing\n", 
					plotFilePath );
				plotfile = stdout;
			}
		}
		else
		{
			plotfile = stdout;
		}
		initplotter();
		numlines = 1;
		drawit( fontname, &xoffset, &yoffset, root );
    	finishplotter();
    	fclose( plotfile );
	}
	else
		syslog(LOG_ERR,"Drawtree: invalid argument list\n");
	return 0;
}
